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Abstract 

We  investigate  the  problem  of  inferring  elastic  moduli  of  nonlinearly  elastic  mem¬ 
branes  from  interior  measurements  of  deformation  and  pressure.  We  begin  by  formu¬ 
lating  a  model  of  membrane  deformation  under  a  vertical  force  where  the  geometry 
of  the  membrane  is  star-like.  The  model  makes  no  specification  of  the  constitutive 
law  by  which  stresses  are  calculated  from  applied  strains.  Under  appropriate  choice 
of  the  constitutive  law  and  simplification  of  the  geometry,  we  show  that  membranes 
of  regular  structure  may  be  homogenized  to  an  axisymmetric  case.  We  then  inves¬ 
tigate  numerical  methods  for  the  resolution  of  the  axisymmetric  model  in  terms  of 
radial  and  vertical  displacement.  Examples  are  given  for  various  boundary  conditions 
and  choices  for  elastic  moduli.  We  present  a  method  by  which  the  moduli  may  be 
accurately  recovered  by  algebraic  calculation  from  knowledge  of  the  displacements  on 
the  boundary  and  interior  of  the  membrane  together  with  measurement  of  the  radial 
stress  at  one  of  the  boundaries. 


RICE  UNIVERSITY 


Recovering  Elastic  Moduli  of  Nonlinear 
Anisotropic  Annular  Membranes  from  Interior 
Measurements  of  Deformation  and  Pressure 

by 

Arthur  David  Cummings 

A  Thesis  Submitted 
IN  Partial  Fulfillment  of  the 
Requirements  for  the  Degree 

Master  of  Arts 


Approved,  Thesis  Committee: 


Associate  Professor  of  Computational  and 
Applied  Mathematics 


Professor  of  Computational  and  Applied 


Professor  of  Computational  and  Applied 
Mathematics  A 


Aladin  M.  Boriek 

Assistant  Professor  of  Medicine 

Baylor  College  of  Medicine 


Houston,  Texas 
December,  1995 


Recovering  Elastic  Moduli  of  Nonlinear 
Anisotropic  Annular  Membranes  from  Interior 
Measurements  of  Deformation  and  Pressure 

Arthur  David  Cummings 


Abstract 

We  investigate  the  problem  of  inferring  elastic  moduli  of  nonlinearly  elastic  mem¬ 
branes  from  interior  measurements  of  deformation  and  pressure.  We  begin  by  formu¬ 
lating  a  model  of  membrane  deformation  under  a  vertical  force  where  the  geometry 
of  the  membrane  is  star-like.  The  model  makes  no  specification  of  the  constitutive 
law  by  which  stresses  are  calculated  from  applied  strains.  Under  appropriate  choice 
of  the  constitutive  law  and  simplification  of  the  geometry,  we  show  that  membranes 
of  regular  structure  may  be  homogenized  to  an  axisymmetric  case.  We  then  inves¬ 
tigate  numerical  methods  for  the  resolution  of  the  axisymmetric  model  in  terms  of 

radial  and  vertical  displacement.  Examples  are  given  for  various  boundary  conditions 

/ 

and  choices  for  elastic  moduli.  We  present  a  method  by  which  the  moduli  may  be 
accurately  recovered  by  algebraic  calculation  from  knowledge  of  the  displacements  on 
the  boundary  and  interior  of  the  membrane  together  with  measurement  of  the  radial 
stress  at  one  of  the  boundaries. 


Acknowledgments 


The  author  would  like  to  take  this  opportunity  to  thank  those  who  have  given  help 
and  support  in  the  completion  of  this  thesis.  First  of  all,  sincere  thanks  to  Steve  Cox, 
my  advisor,  for  the  guidance,  encouragement,  and  critique  he  gave  to  make  this  thesis 
a  better  product.  His  energy  and  insight  formed  the  foundation  of  much  of  my  work. 

Thanks  also  go  out  to  those  who  helped  me  come  to  Rice.  To  Captain  Kevin 
Yeomans  and  Major  Glenn  James,  I  express  my  appreciation  for  their  excitement 
which  kindled  my  interest  in  graduate  study  in  mathematics,  and  for  their  encour¬ 
agement  while  I  sought  to  attend  graduate  school  immediately  after  graduating  from 
the  Air  Force  Academy.  My  gratitude  goes  out  to  Fern  Kinion,  who  worked  a  miracle 
in  getting  the  Air  Force  to  allow  me  to  spend  a  year  and  a  half  here,  and  to  Bob 
Bixby,  who  kept  the  offer  open  as  we  waited  for  the  Air  Force  to  decide. 

I  want  to  express  my  appreciation  for  my  committee  members,  Bill  Symes,  Danny 
Sorensen  and  Aladin  Boriek,  for  their  time  and  energy  in  helping  bring  this  thesis  to 
completion.  Special  acknowledgment  goes  to  Dr.  Boriek  for  the  efforts  he  has  made 
in  acquainting  people  in  this  department  with  his  research  and  exciting  them  to  work 
on  it. 

Finally,  I  want  to  thank  my  wife,  Donna,  for  the  support  and  encouragement  she 
has  given  me  throughout  our  marriage,  and  express  appreciation  for  my  daughter, 
Kayla,  who  is  the  greatest  thing  that  I  have  gained  during  my  time  at  Rice. 


Contents 


Abstract  ii 

Acknowledgments  iii 

List  of  Illustrations  vi 

List  of  Tables  viii 

1  Introduction  1 

2  Models  of  Membrane  Deformation  4 

2.1  General  Membrane  Models .  4 

2.2  A  Model  of  a  Star-like  Membrane  Under  a  Vertical  Force .  5 

2.2.1  Derivation  of  the  Model .  6 

2.2.2  Form  of  the  Constitutive  Law .  15 

2.3  Homogenization  of  the  Model . . .  19 

3  Resolution  of  the  Forward  Problem  22 

3.1  Solving  the  System  of  Equations  for  u"  and  w” .  22 

3.2  Numerical  Solution .  24 

3.2.1  Nonlinear  Shooting  Method .  25 

3.2.2  Nonlinear  Finite  Difference  Method .  25 

3.3  Examples  of  Solutions  to  the  Forward  Problem .  27 

4  A  Direct  Approach  to  the  Inverse  Problem  32 

4.1  Solution  of  the  System  of  Nonlinear  Elasticity  in  Terms  of  the  Moduli  32 

4.1.1  Linear  Isotropic  Constitutive  Law .  33 

4.1.2  Linear  Specially  Orthotropic  Constitutive  Law .  35 

4.2  Examples  of  Recovered  Moduli  .  36 


5  Concluding  Remarks 


40 


V 


Bibliography 


41 


Illustrations 


1.1  A  picture  of  the  canine  diaphragm.  Notice  the  alternating  structure 
of  the  lighter  collagen  bands  with  the  darker  muscle  bundles  near 

around  the  central  tendon .  3 

2.1  An  element  of  the  deformed  membrane  with  associated  unit  vectors. 

Vectors  are  drawn  on  the  positive  faces  with  which  they  are 
associated.  However,  all  vectors  act  through  the  darkened  point 

{r  +  u,0  +  (f)) .  12 

2.2  Representative  volume  element  of  a  unidirectional  fiber  composite 
used  by  Gibson,  as  well  as  Tsai  and  Hahn  in  deriving  expressions  for 

the  effective  moduli  of  the  composite .  17 

3.1  Results  from  solving  the  forward  problem  using  the  information  in 
Table  3.1.  The  solid  line  represents  the  deformed  state,  while  the 

dashed  line  represents  the  flat  reference  state .  28 

3.2  Results  from  solving  the  forward  problem  using  the  information  in 
Table  3.2.  The  solid  line  represents  the  deformed  state,  while  the 

dashed  line  represents  the  flat  reference  state .  30 

3.3  Results  from  solving  the  forward  problem  using  the  information  in 
Table  3.3.  The  solid  line  represents  the  deformed  state,  while  the 

dashed  line  represents  the  flat  reference  state .  31 

3.4  Results  from  solving  the  forward  problem  using  the  information  in 
Table  3.4.  The  solid  line  represents  the  deformed  state,  while  the 

dashed  line  represents  the  flat  reference  state .  31 

4.1  Results  from  recovering  the  elastic  moduli  for  Example  1.  The  solid 
line  represents  the  true  values  of  the  moduli  while  the  stars  represent 
the  recovered  values  of  the  moduli .  37 


vii 

4.2  Results  from  recovering  the  elastic  moduli  for  Example  2.  The  solid 
line  represents  the  true  values  of  the  moduli  while  the  stars  represent 

the  recovered  values  of  the  moduli .  38 

4.3  Results  from  recovering  the  elastic  moduli  for  Example  3.  The  solid 
line  represents  the  true  values  of  the  moduli  while  the  stars  represent 

the  recovered  values  of  the  moduli .  38 

4.4  Results  from  recovering  the  elastic  moduli  for  Example  4.  The  solid 
line  represents  the  true  values  of  the  moduli  while  the  stars  represent 

the  recovered  values  of  the  moduli .  39 


Tables 


2.1  Factors  determining  the  force  vectors  acting  on  the  membrane  ....  13 

3.1  Functions  and  boundary  conditions  characterizing  Example  1 .  28 

3.2  Functions  and  boundary  conditions  characterizing  Example  2 .  29 

3.3  Functions  and  boundary  conditions  characterizing  Example  3 .  29 

3.4  Functions  and  boundary  conditions  characterizing  Example  4 .  30 


1 


Chapter  1 


Introduction 


The  task  of  parameter  identification  from  some  measurement  of  physical  phenomena 
has  application  in  various  fields.  Across  these  fields,  there  seem  to  be  two  major 
approaches  in  measuring  the  physical  occurrences.  Both  approaches  can  be  described 
in  terms  of  the  following  problem: 


-  V-(aVu)  =  /inn  (1.1) 

u  =  g  on  do,  (1-2) 

where  fl  is  a  region  in  9?^.  The  function  a  is  the  parameter  that  we  want  to  identify, 
while  u  is  the  measurable  phenomenon. 

The  first  approach  seeks  to  gain  information  about  a  over  the  entire  domain  by 
measuring  u  only  on  the  boundary  of  the  object.  In  this  case,  the  function  /  would 
be  known  over  the  entire  domain,  and  g  would  be  known  over  the  entire  boundary. 
Such  an  approach  has  appeal  in  that  it  is  not  intrusive  to  the  interior  of  the  object. 
In  application,  such  an  intrusion  may  disturb  the  normal  functioning  of  the  object, 
or  it  may  be  infeasible  for  other  reasons  to  observe  phenomena  on  the  interior  of  the 
object.  Nakamura  and  Uhlmann  [17]  discuss  identification  of  the  Lame  parameters 
which  determine  the  elastic  properties  of  a  linear,  non-homogeneous,  isotropic  elastic 
medium  from  boundary  measurements.  The  success  of  their  method  relies  on  the 
determination  of  a  mapping  A  from  the  Dirichlet  boundary  conditions  of  (1.2)  to 
Neumann  boundary  conditions  of  the  form 

“I  = 

where  n  is  the  outward  unit  normal  to  dD,,  and  a  depends  on  the  Lame  parameters. 
Their  discussion  is  representative  of  this  first  approach  to  parameter  identification. 

A  second  approach  involves  measurement  of  the  phenomena  u  on  the  interior  of  the 
body.  As  in  the  first  approach,  /  and  g  are  known.  However,  here  u  is  known  on  the 
interior  of  H  as  well.  Hence  a  is  the  only  quantity  which  is  unknown.  An  example  of 
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this  approach  is  provided  by  Kohn  and  Lowe  [16]  who  present  a  variational  method 
for  parameter  identification  for  a  similar  problem  in  water  flow.  They  also  review 
other  methods  of  parameter  estimation  which  use  interior  measurements  as  part  of 
their  approach. 

In  this  thesis,  we  will  concern  ourselves  with  inferring  the  elastic  moduli  of  a  non- 
linearly  elastic  membrane.  Our  interest  in  this  application  arises  from  an  effort  to 
better  understand  the  mechanics  of  the  diaphragm  muscle.  It  is  intended  that  this 
work  will  be  integrated  with  ongoing  research  in  in  vivo  animal  physiologic  experi¬ 
mentation  taking  place  in  the  Pulmonary  Section  of  the  Baylor  College  of  Medicine. 
This  experimentation  involves  measurement  of  the  deformation  of  the  diaphragm  un¬ 
der  a  known  pressure.  Knowledge  of  this  deformation  is  analogous  to  knowing  u  on 
both  the  interior  and  the  boundary  of  in  our  example  problem.  Knowledge  of  the 
pressure  is  similar  to  knowing  /.  Boriek,  Wilson,  and  Rodarte  [2]  present  the  details 
of  how  these  measurements  are  made.  Since  their  methodology  does  yield  interior 
measurements  of  deformation,  we  will  focus  on  the  second  approach  to  parameter 
identification. 

An  important  preparation  for  recovering  the  elastic  moduli  is  to  determine  the 
form  of  the  constitutive  law  that  predicts  the  state  of  stress  from  knowledge  of  strain. 
In  the  case  of  the  diaphragm,  this  can  be  extremely  difficult.  Figure  1.1  shows  a 
picture  of  a  canine  diaphragm.  From  this  photograph,  it  is  clear  that  the  diaphragm 
has  a  complicated  structure.  Near  the  chest  wall,  we  can  see  that  radial  bands 
of  collagen  alternate  with  muscle  bundles  as  we  move  around  the  central  tendon. 
Although  the  diaphragm  is  not  homogeneous,  from  this  viewpoint  it  does  seem  to 
have  some  regular  structure.  In  this  thesis,  we  will  develop  a  model  which  will  allow 
us  to  consider  some  aspects  of  this  structure.  We  also  present  a  method  by  which  we 
are  able  to  recover  the  elastic  moduli  of  nonlinear  axisymmetric  membranes  from  the 
interior  measurements  of  the  deformation. 


Figure  1.1  A  picture  of  the  canine  diaphragm.  Notice  the  alternating 
structure  of  the  lighter  collagen  bands  with  the  darker  muscle  bundles  near 

around  the  central  tendon. 
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Chapter  2 

Models  of  Membrane  Deformation 

2.1  General  Membrane  Models 

The  first  step  in  recovering  the  elastic  moduli  which  describe  the  material  properties 
of  the  diaphragm  muscle  is  to  develop  an  appropriate  mathematical  model  which 
describes  the  deformation  of  the  muscle  when  stress  is  applied.  We  have  chosen  to 
model  the  diaphragm  by  an  engineering  structure  known  as  a  membrane.  Since  we 
do  not  wish  to  include  in  our  model  the  strong  central  tendon  of  the  diaphragm,  we 
can  represent  the  diaphragm  by  a  star-like  geometry  with  a  star-like  hole.  We  then 
would  choose  boundary  conditions  at  the  interior  boundary  in  attempt  to  mimic  the 
interface  between  the  diaphragm  and  the  central  tendon.  Further  simplification  to  a 
regular  annulus  would  allow  us  to  consider  models  of  axisymmetric  membranes. 

The  concept  of  a  membrane  is  a  simplification  of  a  more  general  three-dimensional 
engineering  structure  known  as  a  shell.  A  shell  may  be  subjected  to  not  only  stresses, 
but  moments  as  well.  As  a  result,  the  equations  for  the  bending  of  shells  can  be  quite 
difficult  to  formulate.  Green  and  Adkins  [9]  note  that  if  the  changes  in  dimension  vary 
only  slightly  throughout  the  thickness  of  the  deformed  shell,  then  shearing  stresses  in 
this  direction  and  couples  are  then  negligible.  This  is  the  motivation  for  the  theory 
of  elastic  membranes.  Thus  membrane  theory  neglects  shears  normal  to  the  plane  of 
the  membrane,  as  well  as  both  bending  and  twisting  moments.  These  assumptions 
are  reasonable  for  the  passive  diaphragm  which  we  are  trying  to  model.  They  would 
not  hold  for  an  active  diaphragm. 

Green  and  Adkins  formulate  the  theory  of  membranes  independently  of  the  theory 
of  shells.  Antman  [1],  however,  first  derives  equations  for  shells,  then  applies  the 
simplifying  assumptions  which  bring  us  to  membrane  theory.  In  both  cases,  a  great 
deal  of  attention  is  given  to  the  axisymmetric  membrane.  Antman’s  approach  provides 
a  great  deal  of  flexibility  in  specifying  the  undeformed  geometry  of  the  axisymmetric 
membrane,  as  well  as  versatility  in  specifying  the  form  of  the  body  force  to  which  the 
membrane  is  subjected.  This  allows  him  to  consider  not  only  plate  deformations,  but 
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the  deformation  of  spherical  membranes.  Dickey  [6]  presents  in  his  paper  a  model 
for  the  deformation  of  an  initially  flat,  hyperelastic  circular  membrane  under  vertical 
pressure.  His  equations  are  an  instance  of  the  more  general  equations  derived  by 
Antman. 

In  section  2.2  we  develop  a  general  model  for  deformation  of  a  non-homogeneous 
star-like  membrane  subjected  to  a  vertical  pressure  by  generalizing  the  derivation  of 
Dickey’s  model  to  this  domain.  We  do  this  by  first  determining  the  strains  associated 
with  the  deformation.  From  the  resulting  balance  of  forces  in  the  deformed  state,  we 
then  derive  the  equations  of  nonlinear  elasticity.  This  derivation  places  no  restriction 
on  the  form  of  the  stress-strain  relationship.  We  therefore  consider  some  relevant 
forms  of  the  constitutive  law.  In  particular,  we  discuss  a  stress-strain  relationship 
which  allows  periodic  structures  such  as  that  of  the  diaphragm  shown  in  Figure  1.1 
to  be  considered  macroscopically  homogeneous.  In  section  2.3  we  apply  this  consti¬ 
tutive  law  to  homogenize  the  stiffness  of  this  membrane.  With  additional  geometric 
simplifications,  we  then  arrive  at  an  axisymmetric  form  of  the  equations  of  nonlinear 
elasticity. 

2.2  A  Model  of  a  Star-like  Membrane  Under  a  Vertical 
Force  , 

In  this  section,  we  will  develop  a  model  for  a  membrane  which  can  be  expressed 
in  polar  coordinates,  but  is  not  axisymmetric.  Therefore,  the  model  must  allow  for 
the  membrane  to  vary  with  6.  In  doing  so,  we  will  generalize  the  approach  which 
Dickey  [6]  sets  forth  in  the  derivation  of  his  axisymmetric  model. 

We  define  the  reference  state  of  the  membrane  of  thickness  h  to  be,  in  Cartesian, 
coordinates, 

(rcos0,  rsin^,  z{r,d)) 
for  r  and  0  belonging  to  the  domain 

Cl  =  {(r,^)  :  fi{0)  <r  <  f2{0),0  e  [0,27r)}  ^ 

where  fi  and  /2  are  continuous  2%  periodic  functions  of  9.  Note  that  the  reference 
state  need  not  be  flat.  The  vertical  reference  configuration  need  only  be  representable 
as  a  differentiable  function  of  r  and  6,  Due  to  the  application  of  a  vertical  pressure 
of  magnitude  P,  the  membrane  is  deformed  to  the  point 

((r  +  u(r,  9))  cos{9  +  (l){r,  9)),  (r  +  u(r,  9))  sin((9  +  (j){r,  9)),  z{r,  9)  +  w{r,  9)) 
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where  u,  (f>,  and  w  indicate  radial,  angular,  and  vertical  displacements,  respectively. 
We  specify  the  following  displacements  as  boundary  conditions: 

u{h{e),e)  =  a4e),  u{Me),e)  =  ^40), 

t»(/,(9),«)  =  a„(«).  w(h(e),e)=MD),  (2-1) 

2it). 

We  should  note  that  we  may  specify  conditions  on  the  radial  stress  at  the  inner 
and  outer  radius  instead  of  the  radial  displacement.  Introducing  the  function  (Jj. 
representing  radial  stress,  we  can  then  state  the  boundary  conditions  as 

ar{h{d),e)  =  a,(0),  ar{f2{e),e)  =  ^,(9), 

wiMe),e)  =  aM,  Hh{e),e)  =  (2.2) 

(j){r,0)  =  <f>{r, 2Tr). 

To  simplify  notation,  we  will  not  generally  include  the  arguments  of  the  functions 
u,  w,  z,  and  (p,  but  will  assume  them  to  be  understood.  We  will  indicate  the  partial 
derivative  of  one  of  these  functions  with  respect  to  one  of  the  arguments  by  the  func¬ 
tion  subscripted  by  a  comma  and  the  argument  (i.e.  u^r  =  |^).  The  subscripts  r  or  0 
without  a  comma  indicate  direction,  and  are  not  to  be  interpreted  as  differentiation. 


2.2.1  Derivation  of  the  Model 
Derivation  of  Strain  Definitions 


The  first  step  in  the  derivation  of  the  model  is  to  determine  the  equations  which  define 
strain  in  the  radial  and  circumferential  directions.  These  equations  are  derived  from 
differential  considerations  of  the  deformation  function  in  the  appropriate  direction 
using  the  basic  definition  of  strain: 

(deformed  length)  —  (undeformed  length) 
undeformed  length 


strain  = 


Tensile  Strain  in  the  Circumferential  Direction  First,  let  us  consider  the 
deformation  of  an  arc  of  length  z'^g  dO  for  a  given  value  of  r.  Note  that  for  a 

given  r,  the  relationship  above  describes  a  set  of  functions  which  are  parameterized 
by  6  where 

X*  =  [r  +  u{r,9)]cos(0  +  (j){r,0)) 

-  y*  =  [r  +  u(r,  ^)]  sin(^ -f  (/>(r,  0)) 
z*  =  z{r,9)  +  w{r,9). 
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Thus  the  length  of  the  deformed  arc  is  ds*  where 


(ds 


*  \2 


dx* 

He 


+ 


dd 


+ 


Now,  we  have 

dx* 

He 

de 

dz* 

He 


u^e  cos(^  +  <f>)  -  {r  +  u)(l  +  sin(0  + 

Ufi  sin(^  +  (l>)  +  {r  +  u)(l  +  <f>^g)  cos{e  +  (j)) 
z,e  + 


Thus, 

=  Hg  +  (n  +  u)^(l  +  +  iz,g  + 

Therefore,  the  initial  differential  arc  length  of  r  de  is  deformed  to 

+  (r  +  u)2(l  +  +  [z^g  +  w^gf  de. 


Thus  we  can  define  the  tensile  strain  in  the  circumferential  direction  as 

^u'^g  +  (r  +  u)2(l  +  ^^gy  +  {z^g  +  W^eY  de  -  ^JHHHg  de 


Sg  = 


^JH+Hgde 


\ 


Hg  +  (r  +  u)^(l  +  <f>^gY  +  {zfi  +  WfiY 


-  1 


(2.3) 

(2.4) 

(2.5) 


(2.6) 


Tensile  Strain  in  the  Radial  Direction  Similarly,  let  us  consider  the  curve  of 
length  dr,  holding  d  to  be  constant.  We  know  that  the  deformed  arc  length 

is  defined  similarly  by 


{ds:?  = 


Now,  we  have  that 

dx* 

dr 

dy* 

dr 

dz* 

dr 


=  (1  +  u^r)  cos(e  +  ^)  —  (r  +  u)(i)^r  sin(d  +  (j)) 
=  (1  +  U^r)  sin(0  +  </>)  +  (r  +  u)4>,r  cos{e  +  (f)) 
=  z^r  + 


(2.7) 

(2.8) 

(2.9) 
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Thus 


{ds*Y  =  (1  +  U^rY  +  (r  +  +  {z^r  + 

We  also  define  the  tensile  strain  in  the  radial  direction  as 


Er  = 


^J{l  +  U^rY  +  (^  +  ^Y^^r  +  {^,r  +  ^,rY  dv  -  -^1+^ 


1  dr 


yr+^dr 


(1  +  U^rY  +  (^  +  uY(l>'^r  +  {^,r  +  '^,rY 

1  + 


-  1 


(2.10) 


Shear  Strain  We  now  introduce  an  appropriate  set  of  orthogonal  unit  vectors  for 
use  in  describing  the  deformed  state.  These  are 

r  =  cos(0  +  (;^)  z  +  sin(0  +  ^)  j  (2.11) 

0  =  —  sin(0  +  ^)l  +  cos(0  +  (j!>)7  (2.12) 

where  z,J,  and  k  are  the  usual  unit  vectors  in  the  x,y,  and  directions.  Note  that  r 
and  6  are  functions  of  r  and  0. 

With  this  basis  for  the  xy-plane  defined,  we  now  are  interested  in  how  the  deformed 
segments  described  in  sections  2.2.1  and  2.2.1  are  oriented  with  respect  to  these  basis 
vectors.  We  will  consider  only  the  projections  of  these  segments  into  the  icy-plane. 


Shear  Strain  in  the  Circumferential  Direction  The  slope  of  the  segment 
in  section  2.2.1  projected  into  the  a:y-plane  is 
dy*  dy*fd0 

dx*  dx*/d0 

—  sin(^  +  <?i)  +  (r  +  u)(l  +  <f>^e)  cos{0  +  (f)) 

U^g  cos(6>  +  ^)  -  (r  +  'U)(l  +  (j)fi)  sm{0  + 


The  slope  of  0  is 


cos 


(«  +  '!>) 


=  —  coi{0  +  (j>). 


sin(0  +  (f) 

We  will  call  the  angle  between  these  two  slopes  71,  and  will  associate  a  positive  angle 
with  a  deflection  in  the  clockwise  direction.  Thus,  by  use  of  a  trigonometric  identity, 
we  find  that 


tan  7i  = 


—  cot(0  +  (f) 


_  ^ 

dx* 


1  —  cot(0 
u,e 


\i!il 

'  dx* 


(r  +  u)(l  + 


(2.13) 
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Shear  Strain  in  the  Radial  Direction  The  slope  of  the  segment  in  sec¬ 
tion  2.2.1  projected  into  the  a;i/-plane  is 


dx*- 


The  slope  of  r  is 


dy*  I  dr 
dx* (dr 

(1  -1-  u_r)  sin(0  +  <i>)  +  {r  +  u)(j)^r  cos(0  +  (f>) 
(1  +  u^r)  cos{9  +  (j))  —  {r  u)(f>^r  sin(0  +  (f)) 

sm(0  +  4>) 

—  tan(^  T  ^). 


COS(0  +  (f)) 

We  will  call  the  angle  between  these  two  slopes  72,  and  will  associate  a  positive 
angle  with  deflection  in  the  counter-clockwise  direction.  Thus,  by  use  of  the  same 
trigonometric  identity,  we  find  that 

tan(0  -f  4>) 


tan  72  = 


dx* 


1  -}-  tan(^  -f 
{r  -b  u)<i)^r 

1  -j"  r 


(2.14) 


Total  Shear  Strain  The  total  shear  strain  is  the  sum  of  71  and  72.  Thus  we 
may  define  the  shear  strain  as 

u,9  \ 
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arctan 


+  arctan 


1  +  u  - 


(r  -I-  u)(l  -I-  ^,0) 

Applying  the  same  trigonometric  identity,  we  may  rewrite  the  above  equation  as 

Ufi(\  +  u^r)  +  +  ^)^(1  + 


tan  7 


(r  +  u)[(l  + 


(2,15) 


Equilibrium  Conditions  and  Resulting  Equations 

The  deformed  state  of  the  membrane  occurs  when  all  forces  are  in  equilibrium,  i.e., 
when  the  membrane  stresses  balance  the  applied  vertical  pressure.  Having  determined 
the  tensile  and  shear  strains,  we  are  able  to  express  both  the  orientation  and  the  area 
of  the  faces  on  which  the  forces  act.  This  will  enable  us  to  determine  the  equations 
that  must  be  satisfied  in  order  for  equilibrium  to  occur. 


Unit  Normal  and  Tangent  Vectors  to  the  Deformed  Faces  In  order  to  con¬ 
sider  a  balance  of  forces,  we  must  first  understand  the  directions  in  which  the  forces 
act.  Also  we  find  it  useful  to  express  these  vectors  in  terms  of  the  unit  vectors  r  and 
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Unit  Tangent  Vectors  From  equations  (2.3)-(2.5),  we  can  obtain  a  convenient 

representation  for  a  vector  tangent  to  the  face  r  =  constant  as  follows. 

/ 

dx*  _  dy*  _  dz*  -r- 

We  find  the  components  of  in  terms  of  f ,  6,  and  k  to  be 

h  -r  —  Ufi 

ii-0  =  (r  +  m)(1  + 

ti  ■  k  =  z^e  +  w,9 

Thus  the  unit  tangent  vector  to  the  face  r  =  constant  is 

—  +  (r  +  u)(l  +  ^,0)0  +  (^5  + 

J  1  — 


+  (r  +  u)2(l  +  +  {zfi  +  w^sY 


(2.16) 


Similarly,  we  can  define  a  tangent  to  the  face  0  =  constant  using  equations  (2.7)- 
(2.9). 

^2 


dx*  _  dy*  _  dz*  y 


The  components  of  this  vector  in  the  r,  0,  and  k  directions  are 

t2-f  =  1  +  u^r 

h-O  =  {r  +  u)4>^r 

1,2  ■  k  =  +  ^,r 

Thus  the  unit  tangent  vector  to  the  face  6  =  constant  is 

(1  +  U^r)  r  +  (r  +  u)(l>^r  d  +  {z^r  +  to^r)  k 


T2  = 


y(l  +  +  (r  +  uf4>\  +  {2,r  + 


(2.17) 


Unit  Normal  Vectors  We  derive  the  unit  normal  vectors  by  first  noticing  that 
in  the  absence  of  shear  strain,  Ti  is  tangent  to  the  face  r  =  constant  and  normal 
to  the  face  9  =  constant.  Likewise,  r2  is  normal  to  the  face  r  =  constant.  It  is  in 
the  presence  of  shear  strain  that  these  vectors  no  longer  can  double  as  unit  normals. 
Thus  we  shall  derive  the  unit  normals  as  rotations  of  the  unit  tangents  through  the 
angle  7. 
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From  equation  (2.15)  we  can  determine  the  following  relations 
_ +  (r  +  u)(l  + 


sm7  = 


cos  7 


U%{l+U,rr 


(r+uj’a  ’  +  (n  +  u)2(l  d-  (i),ey(t>]r  +  (1  +  4>,eY{^  +  «,r-)^  +  U^e4>‘^ 

_ (1  +  ^,e)(l  +  U,r)  -  U,B(l>,r _ 

’^(r+u)2  ^  +  (?”  +  +  4>,ef‘4''^,T  +  (1  +  4>,eY{^  + 


(2.18) 


(2.19) 

We  obtain  the  unit  normal  vector  to  the  face  r  =  constant  by  rotating  T2  through 
the  angle  —7  in  the  r0-plane. 


/ 


A^i 


cos  7  sin  7 


0\ 


—  sin  7  cos  7  0 

0  01/ 


T, 


(2.20) 


Similarly,  we  may  obtain  the  unit  normal  vector  to  the  face  0  —  constant  by 
rotating  Ti  through  the  angle  7  in  the  r 0-plane. 


/ 


N.= 


cos  7  —  sin  7 


0  \ 


sin  7  cos  7  0 

VO  01/ 


T, 


(2.21) 


It  can  be  shown  that  the  projections  of  A^i  and  Tj  into  the  r0-plane  are  orthogonal, 
as  we  would  expect.  The  same  is  true  of  N2  and  T2- 


Forces  to  Be  Considered  Consider  an  element  of  the  deformed  membrane  as 
shown  in  Figure  2.1.  The  vectors  indicated  are  drawn  on  the  positive  faces.  Table  2.1 
outlines  the  elements  which  determine  the  magnitude  of  the  forces.  Note  that  the 
vertical  pressure  acts  upon  the  undeformed  surface  area  of  the  element,  whereas  the 
stresses  <7^  and  <t$  act  upon  the  deformed  surface  areas  of  their  respective  faces.  Also, 
we  introduce  here  the  shear  stress,  r,  which  is  a  function  of  the  shear  strain,  7,  such 
that  T  =  0  when  7  =  0. 

Equilibrium  Conditions  The  deformed  state  occurs  when  all  forces  are  in  equi¬ 
librium.  This  implies  that  the  sum  of  all  forces  described  in  section  2.2.1  is  the  zero 
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Figure  2.1  An  element  of  the  deformed  membrane  with  associated  unit 
vectors.  Vectors  are  drawn  on  the  positive  faces  with  which  they  are 
associated.  However,  all  vectors  act  through  the  darkened  point  {r  -\-  (f>). 


vector.  The  forces  are  evaluated  on  the  four  lateral  faces  as  well  as  the  top  face. 
Hence,  the  equilibrium  condition  is 

[(<T,  Vi  +  r  Ti)/i  +  (r  +  u)2(l  +  4>^sy  +  {z,e  +  w,8y]  (2.22) 

+  [(<T$  N2  +  T  r2)/i  Ary^(l  +  U^rf  +  (?■  +  ^ 

(Ar)^  —  — 

+  P(rArA^+ =  0 

Dividing  through  by  h  Ar  and  taking  the  limit  as  Ar  and  A9  go  to  zero  yields 

—  |(cr^  Ni  +  T  Ti)^u^q  +  (r  +  n)^(l  +  (j)^eY  +  [z^q  +  (2.23) 

+  ^  {(ere  A^2  +  T  T2)\I{^  +  ^,r)^  +  (r  +  k  =  0 

For  convenience,  we  now  define  the  following  scalars: 


+  ('^  +  +  {^,e  + 

(1  +  U^rY  +  {r  +  uY(tYr  +  (Ar  +  ^,r)^ 

(1  +  U^rY  +  {'>'  +  uY<iYr  +  (Ar  +  ^.r)^ 

■“fe  +  (^  +  ^)^(1  +  (i>,eY  +  (a«  +  '^,oY 


(2.24) 

(2.25) 
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Stress 

Area  of  Face 

Direction 

P 

r  ZV  A0  + 

k 

CT  if> 

h  +  (^  +  «)^(1  + 

Ni 

r 

Ti 

(^d 

h/SryJ{\  +  U^rY  +  (r  +  +  [z^r  +  W,rY 

N2 

T 

T2 

Table  2.1  Factors  determining  the  force  vectors  acting  on  the  membrane 


Substituting  in  the  equations  for  the  unit  normal  and  tangent  vectors,  we  then  may 
rewrite  equation  (2.23)  as 
d 


dr 


{cTrdi  [(1  +  u^r)  COS  7  +  (r  +  sin  7]  r} 

d 


(2.26) 


o  d 

-  —  |(J,.di[(l  +  U,r)sin7  -  (r  +  u)(f>^r  cos  7]^}  +  —  |crrdi(Zr  +  W^r) 


dr 

+  —  +  (r  +  «)(!  +  4>,$)d  +  {z,e  +  w^e)  A:]} 

d 

+  —  {crgd2[u,e  cos  7  -  (r  +  u)(l  +  (j)^e)  sin  7]  r} 
d 


dr 


+  ^  {<yed2[u,d  sin  7  +  (r  +  u)(l  +  (t>^e)  cos  7]  ^  {<75^2(2  $  +  w^e) 

Q  jp'p 

+  ^  {t[(1  +  u^r)r  +  [r  +  u)(l>^rO  +  {z,r  +  n^,T-)^]}  +  =  0 


Q0  LV-  ‘  -,r/’  I  V'  I  ^ 

Recall  that  f  and  0  as  defined  in  equations  (2.11)  and  (2.12)  are  functions  of  r  and  6. 
Thus  we  cannot  simply  pull  them  out  of  the  differentiation  in  order  to  sum  forces  in 
each  of  the  component  directions.  We  must  use  the  product  rule  before  we  can  sum 
forces. 

Note  that  f  and  6  and  their  derivatives  satisfy  the  following  relations 

dr 

dr 
dr 

dr 

m 

de 


=  </’,r  ^ 

=  -<t>,Tr 

=  -(1  ■+</>, 5)  r 
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Using  these  relations,  performing  the  differentiation  on  the  vectors  yields  the  equation 


(Trdi[(l  +  U^r)  COS  7  +  (r  +  u)(f)^r  sin  0  (2.27) 

d 

+  —  {(Trdi[{l  +  U  r)  COS7  +  (r  +  u)(j)r  sin7]}  r 

or 

+  ard\[[\  +  u^r)  sin 7  —  (r  +  u)(j)^r  cos  7](;^,r  r 

.  -  —  {crrdi[{\  +  u,r)sin7  -  (r  +  COS7]}  'O 

d  -  -  d  _ 

+  —  {ardi{z^r  +  U),r)}  k  +  T{r  +  u)Ufi(f)^r  ^  ^  ^ 

d  d 

-  r(r  +  u)(l  +  (l)^0)(j)^rr  +  —  {T(r  +  u)(l  +  (f)^e)}  ^  +  ^,0)}  k 

+  aed2[u^e  cos  7  -  (r  +  u)(l  +  cfi^e)  sin  7](1  +  0 

d 

+  -^{(Ted2[u^0cos'y  -  (rd-u)(l  +  sin7]}  r 

-  a0d2[u^e  sin 7  +  (r  +  u)(l  +  (fy^e)  cos7](l  +  (f)^e)r 

+  ^  {(Jsd2[u,$sin7  +  (r  +  u)(l  +  (l>^e)  cos  7]}  6 

+  —  {aed2{z^e  +  10, e)}  I  +  t(1  +  +  (j)^e)9  +  —  {r(l  +  u,^)}  f 

d  -  d  -  Pr-  - 

-  T{r  +  u)(l  +  (l)A<t>,r  ^  ^  ^  A{^,r  +  ^,r)}  k  +  —  k  =  0 

Separating  into  the  components  of  this  equation  yields  our  system  of  three  partial 
differential  equations  for  the  radial,  vertical,  and  angular  displacements  of  our  non- 
axisymmetric  deformation. 


d 

—  {Ordi  [(1  +  u^r)  cos  7  +  (r  +  u)(f)^r  sin  7]  +  TU^e} 


(2.28) 


d 

+  —  {cred2[u^e  cos  7  -  (r  +  u)(l  +  (fy^e)  sin  7]  +  r(l  +  u^r)} 


+  (Trdi  [(1  +  U,r)  sin  7  -  (r  +  u)^,^  cos  -y]^^r  -  T(r  +  n)(l  +  ()>,0)<t>,r 
-  cr0d2[u^0  sin  7  +  (r  +  t()(l  +  (jy^e)  cos  7](1  +  (fy^g)  -  T{r  +  u)(l  +  (fyA<i>,r  =  0 


d 

{ardi[{l  +  u^r)  sin 7  -  (r  +  u)(jy^r  cos 7]  -  T(r  +  u)(l  +  (fy^g)}  (2.29) 

or 

d 

+  —  {(J0d2[n,0sin7  +  (r  +  u)(l  +  ^,s)cos7]  +  r{r  +  u)<fy^r} 

+  (Trd-I  [(1  +  U^r)  cos  7  +  (r  +  u)(jy^r  sin  7]^,^  +  ^(r  +  u)ufi(f)^r 
+  cred2[ufiCos,^  -  (r  +  u)(l  +  ^,0)sin7](l  +  (jy^g)  +  r(l  +  ■u,r)(l  +  <!>,$)  —  0 
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—  {ardi{z^r  +  w,r)  +  +  w;,e)}  (2.30) 

d  Pt 

+  ^  Wed2{zfi  +  w,9)  +  t(1  +  M,r)}  +  —  =  0 

with  boundary  conditions  as  specified  by  (2.1).  Notice  that  the  result  is  a  nonlinear 
system  of  boundary  value  problems.  As  a  result,  the  questions  of  existence  and 
uniqueness  of  solutions  are  not  easily  answered.  Antman  [1,  pgs  460,  469]  provides 
some  discussion  and  references  concerning  existence  and  uniqueness  in  the  context  of 
three  dimensional  elasticity. 

2.2.2  Form  of  the  Constitutive  Law 

Throughout  our  derivation,  we  have  not  discussed  the  form  of  the  constitutive  law. 
In  considering  the  appropriate  relationship,  it  is  appropriate  to  recall  the  structure 
of  the  problem,  as  illustrated  by  Figure  1.1.  Our  model  certainly  allows  for  the  het¬ 
erogeneous  composition  of  the  material.  We  consider  two  approaches  for  accounting 
for  this  structure  in  the  form  of  the  constitutive  law. 

Treating  the  Membrane  as  Macroscopically  Homogeneous 

Although  we  may  not  wish  to  consider  the  heterogeneity  of  the  membrane  directly, 
we  can  use  some  knowledge  of  the  underlying  structure  to  decide  which  constitutive 
relationship  we  should  consider.  The  entire  diaphragm  should  not  be  expected  to 
respond  the  same  to  radial  strain  as  it  would  to  circumferential  strain,  because  the 
periodic  structure  occurs  in  the  circumferential  direction  only  [2].  This  indicates 
that  we  should  consider  the  membrane  to  have  two  orthogonal  directions  of  elastic 
symmetry — the  radial  direction  and  the  circumferential  direction.  Such  a  material 
is  called  orthotropic.  Furthermore,  since  the  directions  of  symmetry  of  this  material 
coincide  with  the  basis  vectors  which  we  use  in  our  model,  we  consider  the  membrane 
to  be  specially  orthotropic.  This  allows  us  to  express  the  stress-strain  relationship  in 
a  much  simpler  form.  If  we  use  a  linear  constitutive  law,  the  stress-strain  relationship 
is 

(CTj.  \  1  (  Er  l^OrEr 

(70  J  1  ”  I^roEo  Eq 

Here  Er  and  VrO  are  the  Young’s  Modulus  and  Poisson  Ratio  associated  with  strain 
applied  in  the  radial  direction.  Similarly,  Eq  and  UQr  are  the  moduli  associated  with 
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circumferential  stress.  The  matrix  in  this  stress-strain  relationship  is  called  the  stiff¬ 
ness  matrix.  At  first  glance,  there  appears  to  be  four  moduli  that  must  be  recovered. 
However,  under  the  additional  assumption  of  the  existence  of  a  strain  energy  density 
function,  we  may  find  a  relationship  between  the  moduli  which  allows  us  to  compute 
one  in  terms  of  the  others.  A  strain  energy  density  function,  W  is  a.  function  such 
that  the  stresses  are  derived  according  to  the  equations 

dW  dW 
~  Wr'  ~  We' 

A  material  for  which  there  exists  such  a  function  is  called  hyperelastic.  It  can  be 
shown  [8]  that  under  this  assumption,  the  matrix  in  equation  (2.31)  must  be  sym¬ 
metric.  Hence  we  find  that 


Ver  =  ~^^r9-  (2.32) 


Thus  we  have  only  three  elastic  moduli  to  recover  if  we  consider  the  membrane  to  be 
homogeneous  but  orthotropic.  Hyperelasticity  and  orthotropy  are  common  assump¬ 
tions  in  the  literature  regarding  stress-strain  laws  for  muscle  and  other  biological 
tissues  [4,  7,  13,  14,  21,  22]. 

It  should  be  noted  that  a  special  case  of  an  orthotropic  material  is  when  the 
material  responds  to  a  strain  in  exactly  the  same  way,  regardless  of  the  orientation 
of  the  material.  The  material  would  then  be  called  isotropic.  As  we  have  mentioned, 
the  assumption  of  isotropy  is  not  realistic  for  our  problem  because  the  materials 
which  compose  our  membrane  are  not  likely  to  respond  identically  to  strains.  For  an 
isotropic  material,  the  stiffness  matrix  is  simplified  to  yield  the  following  stress-strain 
relationship: 


(T  ip 

o-e 


1  !  E  vE 

TW  1  vE  E 


Sr 

Se 


(2.33) 


Here  E  is  the  Young’s  Modulus  and  v  is  the  Poisson  Ratio.  Thus  we  have  only  two 
elastic  moduli  to  recover  for  an  isotropic  material. 


Treating  the  Membrane  as  a  Unidirectional  Fiber  Composite 

Another  approach  to  considering  the  heterogeneity  of  the  membrane  is  to  look  to  the 
mechanics  of  composite  materials.  Our  observation  of  the  structure  of  the  diaphragm 
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in  chapter  1  directs  our  interest  to  the  study  of  unidirectional  fiber  composites.  Under 
this  approach,  the  muscle  bundles  would  be  considered  the  fibers  which  are  embedded 
in  a  collagen  matrix.  The  biggest  dilference  between  this  approach  and  the  use  of 
an  orthotropic  constitutive  law  is  that  here  we  can  directly  consider  not  only  the 
structure  of  the  muscle,  but  also  the  properties  of  the  individual  materials  which 
make  up  the  membrane,  i.e.  the  muscle  bundles  and  the  collagen  matrix.  Once  we 
have  arrived  at  some  representative  volume  element  of  the  composite  sheet,  we  can  use 
volume  averaged  stresses  which  are  related  to  the  volume  averaged  strains  via  some 
effective  moduli  of  an  equivalent  homogeneous  material.  The  process  of  arriving  at 
these  effective  moduli  can  be  quite  difficult,  however. 

Gibson  [8]  outlines  two  general  methods  to  approach  the  problem.  The  first  is  a 
‘‘theory  of  elasticity”  approach  in  which  the  equations  of  elasticity  must  be  satisfied 
everywhere  in  the  model.  No  simplifying  assumptions  are  made  concerning  stress 
and  strain  distributions.  As  a  result,  solutions  of  any  sort  are  difficult  to  obtain. 
However,  several  analytical  and  numerical  solutions  of  this  type  can  be  found  in  the 
literature  [11,  3,  5,  10]. 


Figure  2.2  Representative  volume  element  of  a  unidirectional  fiber 
composite  used  by  Gibson,  as  well  as  Tsai  and  Hahn  in  deriving  expressions 
for  the  effective  moduli  of  the  composite. 


The  second  approach  is  called  the  “mechanics  of  materiahs”  approach.  Several 
simplifying  assumptions  are  made  in  this  approach  so  that  the  details  of  the  stress 
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and  strain  distributions  are  unnecessary.  In  addition,  the  fiber-packing  geometry  is 
often  unspecified.  These  assumptions  are: 

1.  Fiber  orientation  does  not  change  from  that  indicated  in  the  representative 
volume  element 

2.  Dimensions  of  the  fiber  and  matrix  do  not  change  along  the  length  of  the  element 

3.  The  matrix  and  fiber  are  perfectly  bonded 

4.  Poisson  and  shear  strains  between  the  matrix  and  fiber  are  neglected 

5.  For  longitudinal  (i.e.  radial)  stress,  average  strains  in  the  composite,  fiber,  and 
matrix  are  equal 

6.  For  transverse  (i.e.  circumferential)  stress,  stresses  in  the  composite,  fiber,  and 
matrix  are  equal 


The  first  two  assumptions  allow  us  to  express  the  fiber-packing  geometry  exclusively  in 
terms  of  the  volume  fraction  of  the  fiber,  Vf  and  the  volume  fraction  of  the  matrix, 
where  we  assume  that  vj  ~\-  Vm  —  1-  The  remaining  assumptions  greatly  simplify  the 
state  of  stress  of  the  representative  volume  element.  Using  these  assumptions,  Gibson 
derives  expressions  for  the  effective  moduli  of  the  composite  for  a  material  whose 
representative  volume  element  is  shown  in  Figure  2.2,  where  the  fiber  is  orthotropic 
with  isotropic  matrix.  Tsai  and  Hahn  [20]  derive  the  equivalent  expressions  for  the 
case  where  both  the  matrix  and  the  fiber  are  isotropic.  Since  the  equations  presented 
by  Tsai  and  Hahn  are  for  a  special  case  of  those  derived  by  Gibson,  we  present  those 
equations. 


-F'r  —  ^/(-^/)^  d~ 

1  _ 

Ee  {Ef)e  E^ 


(2.34) 

(2.35) 

(2.36) 


Here  the  subscripts  /  and  m  refer  to  the  fiber  and  matrix,  respectively.  Recall  that 
due  to  the  assumption  of  hyperelasticity,  it  is  not  necessary  to  specify  a  separate 
equation  for  vq^  of  the  composite.  Note  that  E^  and  VrS  computed  by  a  rule  of 
mixtures,  while  Eo  is  computed  by  an  inverse  rule  of  mixtures.  Empirical  data  seem 
to  verify  that  equation  (2.34)  is  a  good  model  of  the  properties  of  the  composite 
material.  However,  equation  (2.35)  often  fails  to  match  empirical  data  sufficiently  [8, 
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20].  This  is  due  to  assumption  6.  Due  to  equilibrium  constraints,  this  assumption 
is  necessarily  true  when  the  matrix  and  fiber  blocks  have  equal  areas  normal  to  the 
circumferential  direction.  Unfortunately,  few  composite  materials  have  this  property. 
There  are  methods  of  correcting  this  error,  but  most  require  additional  knowledge 
of  the  fiber  packing  structure.  Hopkins  and  Chamis  [12]  present  a  refined  model  for 
transverse  properties  using  a  square  fiber-packing  array  and  a  method  of  dividing 
the  representative  volume  element  into  subregions.  Tsai  and  Hahn  [20]  attempt  to 
remedy  the  model  by  allowing  a  Poisson  strain  caused  by  transverse  loading.  For  our 
problem,  we  will  assume  the  simplest  case  to  be  true.  Therefore,  we  will  use  equation 
(2.35)  in  computing  the  transverse  Young’s  modulus  for  the  composite. 

2.3  Homogenization  of  the  Model 

Equations  (2.28)-(2.30)  give  us  a  very  general  model  of  membrane  deformation.  It  is 
within  the  scope  of  this  model  to  allow  the  material  properties  of  the  membrane  to  be 
different  at  every  point  (r,  9)  in  the  domain  of  definition.  We,  however,  are  concerned 
with  a  much  more  regular  and  periodic  structure.  As  we  noted  earlier,  we  assume 
the  muscle  and  connective  tissue  in  the  diaphragm  to  occur  in  a  periodic  fashion.  We 
would  like  to  use  this  knowledge  of  a  much  simpler  structure  to  simplify  the  system  of 
partial  differential  equations  we  developed  in  the  previous  section.  Homogenization 
gives  us  a  tool  to  find  some  average  behavior  over  the  range  of  9  values.  Having 
determined  this  behavior,  we  can  then  consider  the  membrane  to  be  homogeneous 
in  the  circumferential  direction.  The  implications  of  this  homogeneity  would  greatly 
simplify  our  equilibrium  equations. 

Parton  and  Kudryavtsev  [19]  present  formal  mathematical  homogenization  tech¬ 
niques  both  for  periodic  structures,  as  we  are  considering  here,  and  for  regions  with 
a  wavy  boundary.  Oleinik,  Shamev,  and  Yosifian  [18]  discuss  similar  methods  for 
several  problems  in  elasticity.  While  these  methods  have  a  solid  theoretical  basis, 
application  of  these  methods  to  a  system  of  partial  differential  equations  could  be 
difficult  and  tedious.  We  therefore  consider  the  use  of  effective  moduli  to  homogenize 
the  constitutive  law  in  9.  If  we  assume  that  the  structure  of  the  muscle  and  collagen 
is  similar  to  that  shown  in  figure  2.2,  then  we  can  use  the  “mechanics  of  materials” 
approach  discussed  in  section  2.2.2  to  compute  the  effective  moduli.  With  these  effec¬ 
tive  moduli,  we  can  then  consider  the  membrane  to  be  macroscopically  homogeneous 
in  the  circumferential  direction.  We  homogenize  the  boundary  by  simply  redefining 
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our  domain.  If  we  let  fi{0)  =  a  and  /2(^)  =  b  for  a  <  b,  our  domain  becomes 

n  =  {(r,  0)  ;  r  G  (a,  b),  6  G  [0, 27r)} 

If  we  also  redefine  the  function  z  to  be  a  function  of  r  alone,  these  homogenizations 
will  have  the  following  effects  on  our  model; 

1.  There  will  be  no  circumferential  displacement.  Thus  ^  =  0. 

2.  The  radial  and  vertical  displacements,  as  well  as  the  radial  and  circumferential 
stresses,  will  no  longer  vary  with  6.  Thus  u  =  u(r),  w  =  w(r),  (7r  —  o'r{r),  and 
(ye  =  (yeir). 

3.  By  equation  (2.15),  these  first  two  facts  imply  that  the  shear  strain,  7,  is  also 
identically  zero.  Thus  r  =  0. 

4.  Since  there  is  no  more  0  dependence,  all  partial  derivatives  with  respect  to  0 
are  zero.  The  partial  derivatives  with  respect  to  r  become  full  derivatives. 

As  a  consequence  of  fact  4,  we  will  replace  the  subscript  notation  of  differentiation 
with  '  =  Using  these  facts,  the  tensile  strains  defined  in  equations  (2.6)  and  (2.10) 
become 

Se  =  -,  (2.37) 


The  scalars  defined  by  equations  (2.24)-(2.25)  become 

d  y  +  ly  1 

^(1  +  u'Y  +  (z'  +  *^2 

Equation  (2.29)  becomes  a  tautology,  while  the  remaining  equilibrium  equations  sim¬ 
plify  to 


(2.38) 


^  J  (Trjr  +  u)(l  +  U')  \ 

d^\^J{l  +  u'y  +  {z'  +  w'y] 

d  f  crr{r u){z' w')  I 
(iy  I  ^(1  +  u'Y  -f  {z'  +  w'y  j  ^ 


(2.39) 

(2.40) 
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with  either  the  displacement  boundary  conditions 

U{a)  =  Q!„, 
w{a)  =  a^, 

or  the  traction  boundary  conditions 

w{a)  = 

where  a^,,  /3^,  and  Pa  are  constants.  Recall  that  cr^  and  <7$  obey  the  or¬ 

thotropic  constitutive  law  of  equation  (2.31)  with  the  elastic  moduli  then  defined  by 
equations  (2.34)-(2.36).  Note  that  if  we  allow  and  cr^  to  be  general  functions  of 
the  strains,  equations  (2.39)-(2.40)  are  very  similar  to  those  presented  by  Dickey  [6]. 
However,  we  have  now  extended  them  to  allow  non-flat  reference  states  where  the 
reference  vertical  displacement  can  be  written  as  a  function  of  r. 


u{b)  =  Pu, 
w{b)  =  Pu„ 


<Tr(b)  -  Pa, 

w{b)  =  Pu,, 


(2.41) 


(2.42) 
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Chapter  3 

Resolution  of  the  Forward  Problem 

When  we  know  the  properties  of  the  material,  the  vertical  pressure  P,  and  the  mem¬ 
brane  thickness,  /i,  the  only  remaining  unknowns  in  (2.39)-(2.40)  are  radial  displace¬ 
ment,  w,  and  the  vertical  displacement,  w.  We  would  now  like  to  consider  how  to 
solve  these  equations  for  u  and  w.  This  problem  is  called  the  forward  problem.  While 
our  goal  is  to  solve  for  the  material  properties  when  the  displacements  are  known, 
i.e.  the  inverse  problem,  knowledge  of  how  to  solve  the  forward  problem  is  essential 
to  solution  of  the  inverse  problem. 

In  section  3.1  we  look  at  how  to  solve  equations  (2.39)-(2.40)  in  terms  of  the 
second  derivatives  of  the  displacements,  which  is  an  essential  step  for  some  numerical 
methods  which  we  consider  in  section  3,2,  We  then  look  at  some  examples  of  numerical 
solutions  of  equations  (2,39)-(2.40)  with  different  boundary  conditions  and  various 
choices  of  Young’s  moduli  and  Poisson  ratio. 

3.1  Solving  the  System  of  Equations  for  u"  and  w" 

Our  first  step  is  to  rewrite  equations  (2.39)-(2.40)  so  as  to  solve  for  u"  and  w".  We 
then  would  have  a  second  order  system  of  differential  equations  of  the  form 

u"  =  f{r,u,w,u',w') 

n  (  f  f\  \  ‘  / 

w  =  g[r^u^w^u%w  ) 

If  we  can  derive  a  system  of  this  form,  then  the  solution  of  the  forward  problem  will 
simply  be  the  solution  of  a  two  point  boundary  value  problem,  as  is  generally  the 
case  for  a  second  order  differential  equation.  Finding  the  equations  /  and  g  is  not  a 
straightforward  task.  The  second  derivatives  of  u  and  w  arise  in  (2.39)-(2.40)  from 
carrying  out  the  differentiation  in  terms  of  r.  It  is  clear  that  differentiating  terms  such 
as  (1  +  u')  result  in  second  derivatives.  However,  they  also  arise  in  the  differentiation 
of  (Tr,  since  ar  also  depends  on  u'  and  w'.  We  therefore  must  be  able  to  separate  this 
dependence  in  order  to  solve  the  equations  for  u"  and  w". 
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We  see  that  depends  on  u'  and  w'  since  it  is  defined  to  be  a  function  of  Sr-  Recall 
that  in  the  case  of  linear  constitutive  laws,  such  as  those  presented  in  section  2.2.2, 
we  can  write 


where  qf  is  the  first  row  of  the  stiffness  matrix.  Hence 


with 

_  (1  +  u')u"  +  {z'  +  +  w") 

^[1  +  (z')2][(l  +  u'Y  +  {z'  +  w'Y] 
z'z"^{l  +  u'Y  +  [z'  +  w'Y 
[1  +  {z'Yf!'^  ■ 


(3.2) 


Since  the  stiffness  matrix  has  no  dependence  on  u'  or  m',  the  first  term  of  the  above 
sum  has  no  dependence  on  u"  and  w".  The  same  applies  to  any  term  involving  Sg. 
From  (3.2),  we  see  that  is  affine  in  u"  and  w".  Thus  we  can  write 


/  //  I  //  I 


where 


=  qi 


0 

T  I  .,y[l  +  (z'P][(l+«'P  +  (i;'+to')2] 
0 


•5  =  <?1 


,T  I 


( 


£e 


+  9i 


{z'-\-w')z^* 


z^z%^^+u^y+{z^+w^ 


V 


u  _ 

r  r  2 


Using  this  fact,  we  can  then  rewrite  equations  (2.39)-(2.40)  in  the  form 


A 


u 


w 


h, 
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where 
A  = 


<7r(r-i-u)-h(r+u)(li-u')su  _  o-r(r-hu)(l-{-u')^  (Tr(r-i-u)(l-i-u^)(z'-l-w')  I  (r+n)(l+^^)guj 

d  d^  d 

ar(r'\-u)(l-\-u')(z*-\-w‘)  ,  (t-\-u)(z^ -\-w')Su  (7r(r-{‘u)-\-{T-\-u)(z^  r\-w')s^  (Tr{r-\-u){z' 

d^  ^  d  d  d? 


6  = 


Pr 


_  j  I  <7r{r+u){l-\-u'){z'Pw^)z'' 

~  d  ^ 

crr{l+u'){z'+u>')~{r-\-u){P-\-w')s—(Tr{r-{-u)P*  ,  ar{r-\-u){P +xu') 

~  .7  I  rP 


(3.3) 

(3.4) 


with 


d  =  \J{1  +  u'y  +  {z'  +  w'Y  . 

If  this  matrix  A  is  invertible,  then  we  can  write  equations  (2.39)-(2.40)  in  the  form 
of  (3.1),  and  our  second  order  system  of  differential  equations  would  be 


u 


w 


(3,5) 


with  boundary  conditions  (2.41)  or  (2.42).  Of  course,  A  is  invertible  if  and  only  if 
det  A  ^  0.  We  can  see  that  if  r  +  u  =  0,  the  determinant  would  equal  zero.  This 
would  only  happen  if  a  point  in  the  membrane  was  deformed  radially  to  the  origin. 
Under  a  vertical  pressure,  such  a  deformation  is  not  likely  to  occur.  Another  more 
interesting  circumstance  which  would  make  A  singular  is  if  cr^  =  0.  In  the  case  of 
the  linear  constitutive  laws,  this  is  only  possible  if  Sr  and  S$  both  equal  zero,  since 
the  stiffness  matrix  is  necessarily  nonzero.  A  simple  example  of  where  this  would 
occur  is  when  u  —  u'  =  w'  =  z'  —  0.  While  this  is  not  likely  to  occur  in  the  actual 
deformed  state,  it  is  important  in  the  numerical  solution  of  this  problem  to  avoid 
this  singularity — especially  as  an  initial  guess.  Another  place  where  discontinuity 
occurs  is  where  ^(1  +  u'Y  +  {z'  +  vanishes.  Other  than  these  singularities,  A 
is  continuous.  Keller  [15]  gives  conditions  which  ensure  that  solutions  to  this  system 
exist  and  are  unique. 


3.2  Numerical  Solution 

Note  that  the  system  of  equations  (3.5)  is  highly  nonlinear  in  u  and  w.  Thus  we  must 
use  numerical  methods  for  the  solution  of  nonlinear  systems  of  ordinary  differential 
equations.  Two  methods  which  can  be  used  for  problems  of  this  type  are  the  shooting 
method  and  the  finite  difference  method. 


25 


3.2.1  Nonlinear  Shooting  Method 

Using  a  standard  change  of  variables',  (3.5)  can  be  transformed  into  a  first  order 
system  of  ordinary  differential  equations.  After  this  change  of  variables,  a  Runge- 
Kutta  method  can  be  used  to  integrate  the  system  from  a  set  of  initial  values  of  u,  tc, 
w',  and  u)'.  However,  we  do  not  know  that  the  resulting  solution  will  meet  boundary 
conditions  at  endpoint.  Recall  the  boundary  conditions  of  equation  (2.41).  The 
shooting  method  starts  the  Runge-Kutta  method  with  initial  conditions  u{a)  =  au^ 
w{a)  =  and  iteratively  selects  values  of  u'(a)  and  u;'(a)  so  that  the  integration 
returns  values  of  u{b)  and  w{b)  that  are  within  a  certain  tolerance  of  (3^  and 
respectively.  The  code  we  have  implemented  uses  a  secant  update  of  these  initial 
derivative  values.  As  a  result,  quite  a  few  iterations  are  required  for  convergence,  but 
the  method  does  consistently  converge.  One  useful  feature  of  the  shooting  method 
is  that  its  accuracy  is  the  same  as  the  accuracy  of  the  integration  method.  Thus, 
if  we  use  a  fourth-order  Runge-Kutta  method,  our  shooting  method  is  fourth-order 
accurate.  It  is  therefore  relatively  simple  to  improve  the  accuracy  of  the  overall 
shooting  method. 

For  more  general  boundary  conditions,  we  follow  a  similar  scheme.  We  search 
for  initial  values  of  n,  tc,  u',  and  which  yield  a  function  which  matches  the  four 
boundary  conditions.  Hence  we  have  a  system  of  four  nonlinear  equations  in  u(a), 
u;(a),  u'(a),  and  w\a).  This  system  can  then  be  solved  numerically  for  the  proper 
initial  values.  This  approach  allows  for  the  boundary  conditions  to  be  nonlinear 
in  u,u;,u',  and  w\  Therefore,  we  must  use  this  scheme  in  the  case  of  the  traction 
boundary  conditions  of  equation  (2.42),  Keller  [15,  pg  50]  gives  conditions  which 
ensure  that  Newton’s  method  will  converge  uniquely  in  the  numerical  solution  of  this 
system. 

3.2.2  Nonlinear  Finite  Difference  Method 

The  finite  difference  method  uses  finite  difference  approximations  to  the  derivatives 
which  appear  in  (3.5),  resulting  in  a  system  of  nonlinear  equations  in  the  values  of  u 
and  w  at  the  mesh  points  of  the  partition 


a  =  To  <  ri  <  •  •  •  <  =  b. 

We  define  Ui  —  u{ri)  and  Wi  ~  w{ri)  to  be  these  values,  with  boundary  conditions 
uo  =  Wo  =  ayj^  =  I3u^  '^n+i  =  specified  by  (2.41),  For  simplicity,  let  us 
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also  assume  a  uniform  mesh,  i.e.  r^+i  —  ri  =  Ar,  for  i  —  1,  2, . . . ,  n.  Recall  that  with 
equation  (3.5)  we  have  defined  the  functions  /  and  g  of  equation  (3.1)  such  that 

f{r,u,w,u',w') 
g{r,u,w,u',w') 

Applying  the  second  order  finite  difference  approximations  to  the  derivatives,  we  have 
the  following  system  of  2n  equations: 

-«„  +  2ui  -  U2  +  (Ar)  - - - , - - - )  =  0, 

,  o  ,  /  A  ^2  X/  U3-U1W3-WJ 

-ui  +  2u2-U3  +  {Ar)  l{r2,U2,W2, — - — , - - )  =  0, 


j  =  A-H. 


..  ...  '^n—2  ^n— 2  ^ 


-Uyi—2  "h  '^n  “1“  (  At)  /(ry^„X  ?  ^n— 1  ?  ^n— 1 5 

-Un-l  +  2Un  -  (3u  +  {Arff{rn,  Un,  W^, 

-ayj  +  2wi  -W2A  {Aryg{n,ui,wi, 

...  I  O...  ...  I 


h 

h 

'^n—\  l^w  '^n—1 

h 

h 

V/2 

W2 

•)  ^1 7  ^  7 

h 

U3  -  Ui 

2 7  ^2  7  7  7 

W3  -  Wi 

1. 

h 


)  = 


0, 

0, 

0, 

0, 


,  /  A  \2  /  ^n— 2  '^n  '^n—2  \ 

Wn  +  (Ar)  g{rn-\,Ur,-i,Wn^i, - - - , - - - ) 


-Wn-2  +  2u;^_i 

-Wn-1  +  2Wn  -  0W  +  (Ar)^^(  ^ n )  > 


h  ’  h 

l^u  ^n— 1  l^w  ^n— 1 

h  ’  h 


0, 

0. 


Keller  gives  conditions  regarding  continuity  and  boundedness  of  the  derivatives  of  / 
and  g  with  respect  to  u,  w,  u',  and  w'  which  ensure  that  this  system  has  a  unique 
solution  [15,  pg  96]. 

The  more  general  boundary  conditions  of  (2.42)  can  be  met  by  including  equations 
which  measure  the  difference  in  meeting  the  boundary  conditions,  i.e.  equations  of 
the  form 


ar{a)  —  acr  =  0, 

CTrib)  -  =  0, 

m(a)  -  =  0, 

w{b)  -  l^u,  -  0. 
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Recall  that  ar{a)  is  a  nonlinear  function  of  u(a),  u;(a),  u'(a),  and  w\a).  Similarly, 
ar{b)  depends  on  the  values  of  the  same  quantities  at  r  =  6.  Thus  we  must  also 
include  uq,  wq^  ^nd  Wn-\-i  with  the  variables  for  which  we  must  solve.  This 

would  give  us  2(n  +  2)  nonlinear  equations  in  2(n  +  2)  unknowns. 

There  are  many  software  packages  which  provide  nonlinear  equation  solvers  which 
can  be  applied  to  these  systems  to  arrive  at  an  approximate  solution.  For  example, 
the  FORTRAN  subroutines  LMDIFl  and  LMDERl,  part  of  the  MINPACK  collection  of 
optimization  subprograms,  are  available  from  netlib.  Both  subroutines  solve  nonlinear 
systems  of  equations.  They  differ  in  their  requirement  of  Jacobian  information  about 
the  system.  MATLAB  has  a  built-in  function  fsolve  which  also  solves  a  nonlinear 
system  of  equations.  It  can  be  used  with  or  without  providing  a  Jacobian  matrix. 

Convergence  with  the  finite  difference  method  compares  favorably  with  the  shoot¬ 
ing  method,  although  convergence  in  the  problems  with  traction  boundary  conditions 
seems  to  be  much  more  sensitive  to  the  initial  guess.  One  advantage  to  this  method 
over  the  shooting  method  is  that  is  easier  to  “see”  how  changes  in  parameters  may 
affect  the  values  of  u  and  w  at  the  mesh  points.  This  would  be  particularly  useful 
in  solving  a  minimization  problem  where  the  above  functions  are  constraints.  Also, 
it  is  possible  to  formulate  the  finite  difference  method  without  solving  the  system 
of  equations  (3.5).  This  would  avoid  problems  with  singularity  of  the  matrix  A.  A 
disadvantage  of  the  finite  difference  method  is  that  it  is  only  second-order  accurate. 
There  are  fourth-order  finite  difference  formulas,  but  they  are  very  awkward  to  imple¬ 
ment.  Thus  for  the  finite  difference  method,  higher  accuracy  is  better  gained  through 
decreasing  the  step  size,  Ar,  rather  than  increasing  the  accuracy  of  the  method  itself. 

3.3  Examples  of  Solutions  to  the  Forward  Problem 

We  present  for  examples  several  cases  involving  different  boundary  conditions  and 
different  functions  for  the  elastic  moduli.  For  all  cases,  the  data  was  generated  by 
the  nonlinear  shooting  method  using  a  fourth-order  method  to  solve  the  initial  value 
problems.  The  plots  presented  represent  cross-sections  of  the  membrane.  The  entire 
membrane  would  be  generated  by  rotating  the  cross-section  about  r  =  0.  We  assume 
for  all  cases  that  the  magnitude  of  the  vertical  pressure  is  jP  =  100,  and  that  the 
membrane  thickness  is  /i  =  1.  Also  we  specify  the  interval  over  which  r  ranges  to  be 
[0.5, 1.5].  For  simplicity,  we  assume  a  flat  reference  state. 
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Example  1  Our  first  case  is  characterized  by  the  information  shown  in  Table  3.1. 
Since  Er  =  Ee,  this  membrane  is  isotropic.  We  see  in  Figure  3.1  that  the  clamped 
boundary  conditions  cause  a  somewhat  parabolic  deformed  state,  as  we  might  expect. 


Boundary  Conditions 

u(a)  =  0  u(b)  —  0 
w(a)  =  0  w(b)  =  0 

Moduli 

Er(r)  =  97 

Eff(r)  =  97 
i^re{r)  =  0.34 

Table  3.1  Functions  and  boundary  conditions  characterizing  Example  1 


Figure  3.1  Results  from  solving  the  forward  problem  using  the 
information  in  Table  3.1.  The  solid  line  represents  the  deformed  state,  while 
the  dashed  line  represents  the  flat  reference  state. 


Example  2  Our  second  example  uses  different  displacement  boundary  conditions, 
as  well  as  linear  functions  for  the  moduli  Again,  this  membrane  is  isotropic. 
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Boundary  Conditions 

u(a)  =  0  u{h)  =  0 
w{a)  =  1  w{h)  —  0 

Moduli 

Er{r)  =  50  +  50r 
Ee\r)  =  50  +  50r 
Vrdif)  =  0.25  +  O.lr 

Table  3.2  Functions  and  boundary  conditions  characterizing  Example  2 


Example  3  For  the  third  example,  we  present  an  orthotropic  membrane  where  the 
moduli  are  constant.  We  have  chosen  the  values  of  Er  and  Eg  so  that  the  ratio  of  Eg  to 
Er  is  approximately  1.58.  This  is  compatible  with  the  findings  of  Gates,  McCammond, 
Zingg,  and  Kunov  [7]  in  their  testing  of  stiffness  properties  in  the  canine  diaphragm. 
Here  we  also  see  traction  boundary  conditions  in  the  form  of  (2.42).  Since  u  is  not 
required  to  be  zero  at  the  endpoints,  we  see  that  the  deformed  membrane  is  stretched 
in  the  radial  direction. 


Boundary  Conditions 

ar{a)  —  67  0-^(6)  =  90 
r(;(a)  =  1  wih)  =  0 

Moduli 

Er{r)  =  100 

Eg{r)  =  158 
r're(r)  =  0.25 

Table  3.3  Functions  and  boundary  conditions  characterizing  Example  3 


Example  4  Again,  we  have  an  orthotropic  membrane  in  this  case,  and  we  have 
returned  to  the  displacement  boundary  conditions.  However,  we  now  have  functions 
for  the  Young’s  moduli  which  are  sinusoidal.  The  Poisson  Ratio  remains  a  constant 
value.  While  the  difference  between  this  example  and  example  2  may  not  be  visible 
in  the  graphs,  we  expect  that  this  will  be  an  interesting  case  for  the  inverse  problem. 
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Figure  3.2  Results  from  solving  the  forward  problem  using  the 
information  in  Table  3.2.  The  solid  line  represents  the  deformed  state,  while 
the  dashed  line  represents  the  flat  reference  state. 


Boundary  Conditions 

u{a)  —  0  u{b)  =  0 
w[a)  =  1  w(b)  =  0 

Moduli 

Er{r)  =  95  —  40  cos  irr 
Eg{r)  =  120  —  50cos7rr 
ure{r)  =  0.25 

Table  3.4  Functions  and  boundary  conditions  characterizing  Example  4 
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Figure  3.3  Results  from  solving  the  forward  problem  using  the 
information  in  Table  3.3.  The  solid  line  represents  the  deformed  state,  while 
the  dashed  line  represents  the  flat  reference  state. 


Figure  3.4  Results  from  solving  the  forward  problem  using  the 
information  in  Table  3.4.  The  solid  line  represents  the  deformed  state,  while 
the  dashed  line  represents  the  flat  reference  state. 
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Chapter  4 


A  Direct  Approach  to  the  Inverse  Problem 


We  now  turn  our  interest  to  the  solution  of  the  inverse  problem.  We  assume  that  we 
are  given  observed  values  of  the  displacements  u  and  w  at  several  values  of  r  G  [a,  b] 
including  the  endpoints.  This  would  be  equivalent  to  knowing  the  values  at  several 
points  in  (a,  6)  together  with  the  boundary  conditions  of  (2.41).  We  also  know  the 
magnitude  of  the  vertical  pressure  P,  the  membrane  thickness  h,  and  the  reference 
vertical  displacement  z{r).  The  goal  of  the  inverse  problem  is  to  find  the  material 
properties  to  account  for  the  observed  deformation  under  the  known  pressure.  In 
section  4.1  we  present  a  way  of  solving  equations  (2.39)-(2.40)  in  terms  of  the  elastic 
moduli.  We  then  show  examples  of  how  well  we  can  recover  these  moduli  using  this 
approach  in  section  4.2. 


4.1  Solution  of  the  System  of  Nonlinear  Elasticity  in  Terms 
of  the  Moduli 


We  first  consider  solving  the  governing  equations  of  the  deformation  (2.39)-(2.40)  in 
terms  of  the  elastic  moduli  themselves.  Recall  that  the  moduli  occur  in  (2.39)-(2.40) 
in  the  form  of  the  constitutive  law  that  defines  and  cr^.  Thus  we  first  consider  a 
way  of  solving  for  these  stresses  in  terms  of  the  displacements.  Integrating  equation 
(2.40)  from  a  to  r,  we  find  the  following  relation: 

arjr  +  u){z' +  w') 

^(1  +  u'Y  +  (P  +  w'Y  2/i 


where 

^  _  (Trir  +  u){z'  +  w')  ^  Pr"^ 

1^(1  +  u'Y  +  (P  +  w'Y 


(4.2) 


If  we  assume  that  we  may  infer  from  the  observed  values  of  u  and  w  the  values  of 
their  derivatives,  nearly  all  of  the  quantities  in  equation  (4.2)  are  known  at  r  =  a. 
The  only  unknown  value  is  <7,. (a)  since  we  do  not  know  the  elastic  moduli  at  r  =  a. 
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The  strains  on  which  depends  are  known,  however.  Hence,  for  given  initial  values 
of  the  elastic  moduli,  (7  is  a  known  constant  for  this  problem.  However,  if  we  have 
traction  boundary  conditions  of  the  form  of  equation  (2.42)  or  if  the  value  of  cr^(a) 
can  be  measured  experimentally,  then  (7  is  a  known  constant  regardless  of  knowledge 
of  the  moduli  at  r  =  a.  It  is  interesting  to  note  that  we  could  have  integrated  from  r 
to  b.  This  would  have  changed  the  value  of  C  to  the  same  quantity  evaluated  at  r  =  6 
instead  of  r  =  a.  Then  C  would  be  a  known  constant  if  the  all  quantities  are  known 
at  r  =  6.  However,  to  be  concise  we  will  discuss  only  the  case  where  the  quantifies 
are  known  at  r  =  a,  and  C  is  specified  by  equation  (4.2). 

From  equation  (4.1)  we  see  that 

arjr  +  n)(l  +  u')  _  \ u' 

y(l  +  u'Y  +  [z'  +  w'Y  5;'  +  loT  2h  j 
If  z'  +  w'  ^  0,  we  can  substitute  this  relationship  into  equation  (2.39),  which  yields 

dr 


1  +u' 


c 


+  w'  \  2h 

Thus  equations  (4.1)  and  (4.3)  give  us  the  following  system  of  equations. 

I  —  I  1  ,7  f  IJ.,,'  ' 

ere 


(4.3) 


2h  ) 


(4.4) 


(7r  -\-w^ 

1  jL  f  IW  (p  _  Eil\\ 

d  dr  \z^+w'  V  2/i  /  J 

Recall  that  d  =  ^(1  +  u')^  +  {z'  +  and  let  us  call  the  right-hand  side  of  this 
equation  y.  Note  that,  given  the  assumption  that  we  know  the  derivatives  of  u  and 
re,  y  depends  exclusively  on  known  constants  and  a  given  value  of  (7r{a).  It  should 
be  noted  that  this  assumption  of  knowing  the  derivatives  is  not  unreasonable,  as  the 
observed  data  may  be  easily  interpolated  by  a  cubic  spline.  The  resulting  interpolant 
would  be  twice  differentiable,  allowing  us  to  construct  the  right-hand  side  of  this 
equation.  Having  arrived  at  this  equation,  we  must  now  consider  the  form  of  the 
constitutive  law  in  order  to  recover  the  elastic  moduli.  We  will  consider  the  linear 
isotropic  and  specially  orthotropic  cases  described  in  Section  2.2.2. 


4.1.1  Linear  Isotropic  Constitutive  Law 

For  an  isotropic  material,  the  linear  constitutive  law  is  as  shown  in  equation  (2.33). 
Using  this  relationship,  equation  (4.4)  becomes 

(Sr  +  lySe) 

^  {Zq  +  vZr) 


l-i/2 
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Thus  we  have  two  equations  and  two  unknowns-the  Young’s  Modulus,  E,  and  the 
Poisson  ratio,  u.  We  expect  there  to  be  a  unique  solution.  In  order  to  see  this,  we 
introduce  the  following  change  of  variables. 


6  = 


6  = 


E 


1  —  Z/2 

vE 

1  —  1/2 


If  E  ^  0,  we  can  change  variables  back  using  the  relationships 


u  = 

E  = 


k 

6 


Since  we  usually  expect  ^  0,  this  is  a  valid  change  of  variables.  Thus  we  can  write 


/  £r  £e 
\  £0  £r 


x  =  y 


(4.5) 


where 


Let  us  call  the  matrix  in  equation  (4.5)  G,  Clearly,  G  is  invertible  whenever  ^  £q  . 
Recalling  equations  (2.37)-(2.38),  we  see  that  this  condition  is  equivalent  to 


7^  r 


\j{  \  +  +  {z'  4.  w'Y 


(4.6) 


Note  that  G  depends  only  on  the  observed  displacements  and  their  derivatives.  Hence, 
if  (4.6)  holds,  the  solution  vector 


x  =  G 


is  determined  uniquely  for  a  known  value  of  ar{a).  Thus  we  can  calculate  algebraically 
the  functions  E{r)  and  u{r)  over  the  entire  interval,  if  we  have  the  necessary  values  at 
r  =  a  and  if  there  is  no  r  G  (a,  b]  such  that  w'{r)  =  ~z\r).  Even  if  there  is  such  an  r, 
we  may  still  be  able  to  determine  the  elastic  moduli  except  in  a  small  neighborhood 
of  this  r. 
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4.1.2  Linear  Specially  Orthotropic  Constitutive  Law 

The  form  of  the  linear  specially  orthotropic  constitutive  relationship  is  given  in  equa¬ 
tion  (2.31).  Using  this  equation,  we  can  write  equation(4.4)  as 


l-^rO^er 

Ef) 


(Sr  +  t'BrSe) 

(Se  +  I'roSr) 


Using  equation  (2.32)  resulting  from  the  assumption  of  hyperelasticity,  we  may  rewrite 
this  equation  as 

/  2  (ErSr  +  VroEoSe)  \ 

I  re  I  —  />/ 


{Sff 


where  we  define 


Thus  there  are  three  elastic  moduli  to  recover  from  this  problem.  However,  since  we 
have  only  two  equations,  we  can  not  expect  to  recover  all  three.  Therefore,  let  us 
assume  that  the  major  Poisson  ratio,  Urs  is  a  known  constant.  Then  we  have  only  the 
longitudinal  and  transverse  Young’s  moduli,  Er  and  E$  to  determine.  We  therefore 
introduce  the  following  change  of  variables: 

Er 

Ee 


6  = 


1  -  P^^re 


If  Er  7^  0,  we  can  change  variables  back  by  using  the  relationships 

6 

Er  =  6  (l  -  P^re) 

Ee  =  6  (l  -  P^le) 

Again,  we  expect  that  Er  0.  Thus  we  consider  this  a  valid  change  of  variables. 
Using  these  new  variables,  we  can  write 


where 


Gx  =  y 

Sr  l^r  9  S$ 

0  VrsSr  +  £$ 
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This  matrix  is  singular  when  =  0  or  when  UroSr  =  —£e>  Otherwise,  we  have  a 
solution  vector 

which  is  uniquely  determined  for  a  given  value  of  and  the  constant  Thus 

if  there  is  no  r  G  such  that  w\t)  =  — 2r'(r),  we  can  calculate  algebraically  the 
functions  Er{r)  and  EQ{r),  Once  these  functions  are  recovered,  we  would  like  to  use 
equations  (2.34)“(2.35)  to  recover  the  moduli  of  the  matrix  and  fiber.  If  the  fiber 
and  matrix  are  both  isotropic  and  the  volume  fractions  vj  and  are  known,  these 
equations  form  a  nonlinear  system  of  equations  in  two  unknowns  which  we  should  be 
able  to  solve  for  Ej  and  E^^ 

4.2  Examples  of  Recovered  Moduli 

We  present  as  examples  results  from  recovering  the  moduli  from  the  four  examples 
in  section  3.3.  In  every  case,  we  use  the  data  for  u  and  w  with  their  first  and  second 
derivatives  which  was  generated  by  the  forward  solution  of  the  problem  using  the 
shooting  method.  We  assume  that  is  known  for  all  cases.  In  example  3,  this 

simply  involves  the  boundary  conditions.  For  other  cases,  we  calculate  this  value 
using  known  initial  values  of  the  elastic  moduli. 

Example  1  The  parameters  of  this  problem  were  specified  in  Table  3.1.  Since  this 
is  an  isotropic  problem,  we  expect  to  recover  both  the  Young’s  modulus  and  the 
Poisson  ratio.  We  can  see  in  Figure  3.1  that  there  is  a  critical  point  where  ~  0. 
Thus  we  expect  to  lose  accuracy  in  recovering  the  moduli  near  this  critical  point. 
Figure  4.1  shows  the  results  of  this  direct  method  using  ar{a)  =  78.0146.  We  clearly 
see  a  spike  in  both  moduli  near  this  critical  point.  However,  the  effect  seems  to  be 
localized,  since  the  recovered  values  seem  to  be  very  close  to  the  true  values  except 
for  only  a  few  points.  Thus,  except  for  a  small  neighborhood  near  the  critical  point 
we  are  able  to  recover  the  moduli  over  the  entire  interval  with  good  accuracy. 

Example  2  In  Figure  3.2  we  can  see  that  we  should  not  have  the  same  difficulty 
with  critical  points.  So  we  hope  to  find  the  correct  values  for  E  and  u  at  all  points 
on  the  interval.  Using  crr{a)  =  44.5845  we  get  the  results  shown  in  Figure  4.2  Clearly, 
all  recovered  values  are  very  close  to  the  true  values. 
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Figure  4.1  Results  from  recovering  the  elastic  moduli  for  Example  1.  The 
solid  line  represents  the  true  values  of  the  moduli  while  the  stars  represent 
the  recovered  values  of  the  moduli. 


Example  3  Since  this  problem  is  an  orthotropic  membrane,  we  can  only  recover 
and  Ee-  We  assume  that  the  value  of  Urs  =  0.25  is  known.  For  this  example,  we  use 
the  value  of  (7r{a)  provided  by  the  boundary  conditions  shown  in  Table  3.3.  Again, 
we  see  very  good  results  in  recovering  the  true  moduli. 

Example  4  For  this  orthotropic  problem,  we  again  assume  that  the  value  of  Uro  = 
0.25  is  known.  Recall  that  the  functions  describing  the  moduli  as  shown  in  Table  3.4 
are  sinusoidal.  Figure  4.4  shows  the  results  of  recovering  these  moduli  using  <7^(0)  = 
49.6505.  There  appears  to  be  slightly  more  error  in  recovering  these  moduli  than  in 
the  other  examples.  This  may  be  due  to  the  wider  variation  in  moduli  values  than 
has  been  present  in  the  other  examples 
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Figure  4.2  Results  from  recovering  the  elastic  moduli  for  Example  2.  The 
solid  line  represents  the  true  values  of  the  moduli  while  the  stars  represent 
the  recovered  values  of  the  moduli. 


Figure  4.3  Results  from  recovering  the  elastic  moduli  for  Example  3.  The 
solid  line  represents  the  true  values  of  the  moduli  while  the  stars  represent 
the  recovered  values  of  the  moduli. 


Circumferential  Young’s  Modulus  Radial  Young’s  Modulus 
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Figure  4.4  Results  from  recovering  the  elastic  moduli  for  Example  4.  The 
solid  line  represents  the  true  values  of  the  moduli  while  the  stars  represent 
the  recovered  values  of  the  moduli. 
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Chapter  5 


Concluding  Remarks 


In  this  thesis,  we  have  focused  our  attention  on  recovering  the  elastic  moduli  of 
axisymmetric,  nonlinearly  elastic  membranes  which  follow  linear  constitutive  laws. 
We  have  shown  that  with  knowledge  of  the  deformation  at  the  boundaries  and  on  the 
interior,  we  can  determine  the  elastic  moduli  quite  accurately,  if  we  also  know  the 
radial  stress  at  either  boundary.  Practical  application  of  this  method  of  recovering 
the  moduli  has  the  following  steps: 

1.  Obtain  measurements  of  vertical  and  radial  displacement  of  the  deformed  mem¬ 
brane,  as  well  as  measurements  of  the  vertical  pressure,  membrane  thickness, 
and  the  radial  stress  at  the  inner  radius. 

2.  Interpolate  the  measurements  of  displacement  with  a  cubic  spline  or  other  in- 
terpolant  which  has  at  least  two  continuous  derivatives. 

3.  At  each  point  where  the  moduli  are  to  be  determined,  solve  equation  (4.4)  for 
the  moduli  using  the  interpolants  to  obtain  the  necessary  derivatives.  This 
solution  may  involve  a  change  of  variables  to  make  (4.4)  a  linear  system. 

Thus  the  determination  of  these  parameters  comes  from  direct  solution  of  algebraic 
equations.  A  possible  area  for  additional  research  with  this  method  would  be  in 
solving  an  optimization  problem  for  the  correct  value  of  radial  stress  at  the  boundary. 
The  solution  to  such  a  problem  might  eliminate  the  necessity  for  knowing  this  value 
prior  to  determining  the  moduli. 

We  have  also  presented  a  model  for  star-like  membrane  deformation  which  is 
much  more  general  than  the  axisymmetric  case.  This  model  is  general  enough  to 
be  applied  to  the  diaphragm  with  far  fewer  simplifying  assumptions.  It  is  hoped 
that  further  research  using  this  model  will  allow  more  accurate  estimation  of  the 
parameters  describing  deformation  of  the  diaphragm. 
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